Compiled under R version 4.2.0 (2022-04-22)
WARNING: edit the working directory to your preferred folder.
This document details all analyses performed in R for the
study:
Legendre, L. J., S. Choi, and J. A. Clarke. The use of inconsistent
terminology for reptile eggshell traits affects the outcome of
evolutionary analyses. Journal of Anatomy 241:641–666.
For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)
tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")
data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match
## character(0)
rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)
thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)
# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))
# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
phytools – 1000 iterations, using
AIC to select best model (code modified from Liam Revell, see
here)x<-thisstudy
aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
d.aic<-aic-min(aic)
exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010427592 0.0005213796 0.0005213796
## Semi-rigid 0.0005213796 -0.0010427592 0.0005213796
## Soft 0.0005213796 0.0005213796 -0.0010427592
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010972800 0.0004587151 0.0006385650
## Semi-rigid 0.0004587151 -0.0007221278 0.0002634128
## Soft 0.0006385650 0.0002634128 -0.0009019778
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.000000e+00 0.0000000000
## Semi-rigid 0.0065561168 -1.015329e-02 0.0035971777
## Soft 0.0004585074 4.433192e-05 -0.0005028393
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
logL
## ER SYM ARD
## -56.60363 -56.31510 -44.90060
AIC<-mapply(aic,logL,c(1,3,6))
AIC
## ER SYM ARD
## 115.2073 118.6302 101.8012
AIC.W<-aic.w(AIC)
AIC.W
## ER SYM ARD
## 0.0012254128 0.0002213112 0.9985532760
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 1 0 999
o1<-make.simmap(tree,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.000000e+00 0.0000000000
## Semi-rigid 0.0065561168 -1.015329e-02 0.0035971777
## Soft 0.0004585074 4.433192e-05 -0.0005028393
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
o2<-make.simmap(tree,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010427592 0.0005213796 0.0005213796
## Semi-rigid 0.0005213796 -0.0010427592 0.0005213796
## Soft 0.0005213796 0.0005213796 -0.0010427592
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
treesstudy<-c(o1,o2)
objstudy<-describe.simmap(treesstudy)
plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 208!
=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.
x<-norell
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.001408430 0.000704215 0.000704215
## Semi-rigid 0.000704215 -0.001408430 0.000704215
## Soft 0.000704215 0.000704215 -0.001408430
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013987686 0.0004876351 0.0009111335
## Semi-rigid 0.0004876351 -0.0010907495 0.0006031144
## Soft 0.0009111335 0.0006031144 -0.0015142479
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.007810285 -0.012947744 0.005137459
## Soft 0.001558250 0.001343177 -0.002901427
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
logL
## ER SYM ARD
## -70.01788 -69.49006 -56.79592
AIC<-mapply(aic,logL,c(1,3,6))
AIC
## ER SYM ARD
## 142.0358 144.9801 125.5918
AIC.W<-aic.w(AIC)
AIC.W
## ER SYM ARD
## 2.685989e-04 6.162366e-05 9.996698e-01
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 0 0 1000
treesnorell<-make.simmap(tree,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.007810285 -0.012947744 0.005137459
## Soft 0.001558250 0.001343177 -0.002901427
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
objnorell<-describe.simmap(treesnorell)
plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Completely different result – let us check which taxa are different from the previous ASR…
datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]
# Only 7 taxa are different – two non-avian dinosaurs and five turtles
# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0014260407 0.0007130204 0.0007130204
## Semi-rigid 0.0007130204 -0.0014260407 0.0007130204
## Soft 0.0007130204 0.0007130204 -0.0014260407
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013980396 0.0005877562 0.0008102834
## Semi-rigid 0.0005877562 -0.0013213236 0.0007335675
## Soft 0.0008102834 0.0007335675 -0.0015438509
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.000000000 0.000000000
## Semi-rigid 0.0076499485 -0.012737081 0.005087133
## Soft 0.0005584858 0.000722445 -0.001280931
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 0 0 1000
treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.000000000 0.000000000
## Semi-rigid 0.0076499485 -0.012737081 0.005087133
## Soft 0.0005584858 0.000722445 -0.001280931
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)
All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.
x<-legendre
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007530010 0.0003765005 0.0003765005
## Semi-rigid 0.0003765005 -0.0007530010 0.0003765005
## Soft 0.0003765005 0.0003765005 -0.0007530010
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006161968 0.0000000000 0.0006161968
## Semi-rigid 0.0000000000 -0.0005646245 0.0005646245
## Soft 0.0006161968 0.0005646245 -0.0011808213
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004310504 0.000000000 0.0004310504
## Semi-rigid 0.0036277349 -0.007172395 0.0035446601
## Soft 0.0003893632 0.000000000 -0.0003893632
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
logL
## ER SYM ARD
## -45.62898 -43.24564 -41.32700
AIC<-mapply(aic,logL,c(1,3,6))
AIC
## ER SYM ARD
## 93.25795 92.49129 94.65400
AIC.W<-aic.w(AIC)
AIC.W
## ER SYM ARD
## 0.3372989 0.4948726 0.1678286
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 337 495 168
l1<-make.simmap(tree,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007530010 0.0003765005 0.0003765005
## Semi-rigid 0.0003765005 -0.0007530010 0.0003765005
## Soft 0.0003765005 0.0003765005 -0.0007530010
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
l2<-make.simmap(tree,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006161968 0.0000000000 0.0006161968
## Semi-rigid 0.0000000000 -0.0005646245 0.0005646245
## Soft 0.0006161968 0.0005646245 -0.0011808213
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
l3<-make.simmap(tree,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004310504 0.000000000 0.0004310504
## Semi-rigid 0.0036277349 -0.007172395 0.0035446601
## Soft 0.0003893632 0.000000000 -0.0003893632
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
treeslegendre<-c(l1,l2,l3)
objlegendre<-describe.simmap(treeslegendre)
plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.
The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.
To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.
remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))
# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008859112 0.0004429556 0.0004429556
## Semi-rigid 0.0004429556 -0.0008859112 0.0004429556
## Soft 0.0004429556 0.0004429556 -0.0008859112
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -7.746782e-04 4.816389e-05 0.0007265143
## Semi-rigid 4.816389e-05 -7.860575e-04 0.0007378936
## Soft 7.265143e-04 7.378936e-04 -0.0014644079
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001530468 0.0000000000 0.0001530468
## Semi-rigid 0.0042378920 -0.0042378920 0.0000000000
## Soft 0.0020047508 0.0008991591 -0.0029039099
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 233 121 646
s1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008859112 0.0004429556 0.0004429556
## Semi-rigid 0.0004429556 -0.0008859112 0.0004429556
## Soft 0.0004429556 0.0004429556 -0.0008859112
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
s2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -7.746782e-04 4.816389e-05 0.0007265143
## Semi-rigid 4.816389e-05 -7.860575e-04 0.0007378936
## Soft 7.265143e-04 7.378936e-04 -0.0014644079
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
s3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001530468 0.0000000000 0.0001530468
## Semi-rigid 0.0042378920 -0.0042378920 0.0000000000
## Soft 0.0020047508 0.0008991591 -0.0029039099
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
treenewsr<-c(s1,s2,s3)
objsr<-describe.simmap(treenewsr)
plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.
x[c(21,22)]<-"Soft"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008789375 0.0004394688 0.0004394688
## Semi-rigid 0.0004394688 -0.0008789375 0.0004394688
## Soft 0.0004394688 0.0004394688 -0.0008789375
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.000754388 0.0000000000 0.0007543880
## Semi-rigid 0.000000000 -0.0005361558 0.0005361558
## Soft 0.000754388 0.0005361558 -0.0012905438
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001518191 0.0000000000 0.0001518191
## Semi-rigid 0.0024414178 -0.0024414178 0.0000000000
## Soft 0.0021283836 0.0005467983 -0.0026751820
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 59 151 790
so1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008789375 0.0004394688 0.0004394688
## Semi-rigid 0.0004394688 -0.0008789375 0.0004394688
## Soft 0.0004394688 0.0004394688 -0.0008789375
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
so2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.000754388 0.0000000000 0.0007543880
## Semi-rigid 0.000000000 -0.0005361558 0.0005361558
## Soft 0.000754388 0.0005361558 -0.0012905438
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
so3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001518191 0.0000000000 0.0001518191
## Semi-rigid 0.0024414178 -0.0024414178 0.0000000000
## Soft 0.0021283836 0.0005467983 -0.0026751820
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
treenews<-c(so1, so2, so3)
objs<-describe.simmap(treenews)
plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Same result as with semi-rigid coding, unsurprisingly.
x[c(21,22)]<-"Hard"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007755266 0.0003877633 0.0003877633
## Semi-rigid 0.0003877633 -0.0007755266 0.0003877633
## Soft 0.0003877633 0.0003877633 -0.0007755266
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006356856 0.000000000 0.0006356856
## Semi-rigid 0.0000000000 -0.000568657 0.0005686570
## Soft 0.0006356856 0.000568657 -0.0012043425
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004481032 0.000000000 0.0004481032
## Semi-rigid 0.0036159226 -0.007150659 0.0035347368
## Soft 0.0003936984 0.000000000 -0.0003936984
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
## ER SYM ARD
## 350 494 156
h1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007755266 0.0003877633 0.0003877633
## Semi-rigid 0.0003877633 -0.0007755266 0.0003877633
## Soft 0.0003877633 0.0003877633 -0.0007755266
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
h2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006356856 0.000000000 0.0006356856
## Semi-rigid 0.0000000000 -0.000568657 0.0005686570
## Soft 0.0006356856 0.000568657 -0.0012043425
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
h3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004481032 0.000000000 0.0004481032
## Semi-rigid 0.0036159226 -0.007150659 0.0035347368
## Soft 0.0003936984 0.000000000 -0.0003936984
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333
## Done.
treenewh<-c(h1, h2, h3)
objh<-describe.simmap(treenewh)
plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.
To check how strongly are these results influenced by branch length
information, we can replicate these analyses using maximum parsimony,
which does not consider branch length information, using
castor.
studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)
MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
# ancestral states
MPS2<-MPS$ancestral_likelihoods
rownames(MPS2)<-c(209:374); colnames(MPS2)<-c("soft","semi-rigid","hard")
norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)
MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)
# ancestral states
MPN2<-MPN$ancestral_likelihoods
rownames(MPN2)<-c(209:374); colnames(MPN2)<-c("soft","semi-rigid","hard")
legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)
MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
# ancestral states
MPL2<-MPL$ancestral_likelihoods
rownames(MPL2)<-c(209:374); colnames(MPL2)<-c("soft","semi-rigid","hard")
xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)
# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")
# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")
# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")
treedi<-multi2di(tree, random=FALSE)
# To visualize node numbers on the tree
ggtree(treedi) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)
contMap assumes a Brownian motion
model)phylosig(treedi, AET, method="lambda", test=TRUE)
##
## Phylogenetic signal lambda : 0.999934
## logL(lambda) : -315.228
## LR(lambda=0) : 257.077
## P-value (based on LR test) : 7.44311e-58
phylosig(treedi, RET, method="lambda", test=TRUE)
##
## Phylogenetic signal lambda : 0.999934
## logL(lambda) : -231.158
## LR(lambda=0) : 207.137
## P-value (based on LR test) : 5.78775e-47
Very strong signal for both absolute and relative eggshell thickness.
phyloplotAET<-contMap(treedi, AET, plot=FALSE)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, AET))
## Ancestral character estimates using fastAnc:
## 209 210 211 212 213 214
## 29.293763 29.293763 26.976449 20.021097 37.072997 57.381022
## 215 216 217 218 219 220
## 59.447111 38.957535 37.945828 62.327061 63.691479 60.852431
## 221 222 223 224 225 226
## 18.697740 17.464853 16.996346 21.391647 14.489654 17.150500
## 227 228 229 230 231 232
## 16.912464 16.433291 14.515531 16.129348 8.218623 8.892389
## 233 234 235 236 237 238
## 17.599849 17.653945 11.371660 11.392832 11.093155 7.684835
## 239 240 241 242 243 244
## 18.184274 27.793329 11.641005 6.394956 5.262023 11.870736
## 245 246 247 248 249 250
## 10.875623 12.261946 12.123157 14.012763 13.895692 14.226702
## 251 252 253 254 255 256
## 16.557906 31.172542 39.050721 131.436576 544.304122 573.403854
## 257 258 259 260 261 262
## 610.641866 607.826064 498.405442 134.897657 140.865662 142.716126
## 263 264 265 266 267 268
## 139.304069 142.860600 61.041297 154.746794 178.179127 198.542856
## 269 270 271 272 273 274
## 196.759526 167.261430 145.626598 217.210263 259.744683 324.035230
## 275 276 277 278 279 280
## 36.490983 252.266971 278.041830 413.790201 623.886942 366.312141
## 281 282 283 284 285 286
## 459.981899 450.032672 468.579008 478.030338 31.529177 2.910210
## 287 288 289 290 291 292
## 0.014539 4.051692 11.279474 3.368817 38.203136 115.277972
## 293 294 295 296 297 298
## 115.277972 1422.220619 1422.220619 115.277972 1810.898354 1810.898354
## 299 300 301 302 303 304
## 1423.663197 1481.546260 29.651406 14.412124 72.463867 1.797803
## 305 306 307 308 309 310
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625
## 311 312 313 314 315 316
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625
## 317 318 319 320 321 322
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625
## 323 324 325 326 327 328
## 319.220006 479.414045 856.300831 856.300831 856.300831 535.750860
## 329 330 331 332 333 334
## 538.224713 538.224713 560.814605 1288.673744 571.363791 610.271423
## 335 336 337 338 339 340
## 610.271423 610.271423 610.271423 1392.053831 1101.924581 544.995844
## 341 342 343 344 345 346
## 544.995844 544.995844 681.611306 681.611306 420.284152 420.284152
## 347 348 349 350 351 352
## 420.284152 401.473506 401.473506 401.473506 401.473506 385.047899
## 353 354 355 356 357 358
## 394.698758 483.723330 483.723330 483.723330 483.723330 483.723330
## 359 360 361 362 363 364
## 488.479400 483.723330 672.731467 688.660804 743.447690 750.028948
## 365 366 367 368 369 370
## 642.663152 1304.111328 319.380290 778.921709 908.582073 820.930413
## 371 372 373 374 375 376
## 423.106033 255.141714 300.722998 343.668052 298.875220 290.389862
## 377 378 379 380 381 382
## 222.173714 217.770862 177.872861 166.751051 161.039115 229.707740
## 383 384 385 386 387 388
## 234.260007 236.237640 233.565993 232.205729 231.164466 196.488111
## 389 390 391 392 393 394
## 216.864203 210.923903 208.771974 210.728366 206.860831 201.759935
## 395 396 397 398 399 400
## 251.229188 289.706525 316.999616 253.504459 257.542381 243.366955
## 401 402 403 404 405 406
## 151.843383 155.880089 164.029302 144.791426 128.410742 106.962875
## 407 408 409 410 411 412
## 100.299131 118.844349 109.671081 73.950228 103.670078 99.287810
## 413 414 415
## 98.491772 103.798544 101.498594
Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.
Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).
phyloplotRET<-contMap(treedi, RET, plot=FALSE)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, RET))
## Ancestral character estimates using fastAnc:
## 209 210 211 212 213 214 215 216
## 3.435293 3.435293 4.047709 7.765196 20.362606 37.938372 39.824787 25.690342
## 217 218 219 220 221 222 223 224
## 39.657428 42.678094 44.856834 53.006727 7.508812 7.363644 7.173682 6.553561
## 225 226 227 228 229 230 231 232
## 9.463168 7.383174 7.655193 9.189493 11.648581 13.022670 12.734649 13.646320
## 233 234 235 236 237 238 239 240
## 16.530771 20.145434 21.874698 24.370253 23.820866 14.769114 1.647011 1.003134
## 241 242 243 244 245 246 247 248
## 1.610621 0.936574 0.860310 1.828622 1.596660 2.012543 3.474795 1.682980
## 249 250 251 252 253 254 255 256
## 1.180146 1.164828 2.085553 3.133934 2.630136 6.622457 10.830493 10.810476
## 257 258 259 260 261 262 263 264
## 10.535636 10.641238 9.256621 6.895434 8.089943 10.710668 10.792255 12.171031
## 265 266 267 268 269 270 271 272
## 2.043677 11.100157 10.747658 10.935592 10.272658 9.993139 8.767617 10.221160
## 273 274 275 276 277 278 279 280
## 7.286643 7.960715 2.157566 4.753714 4.814619 5.477856 6.779844 4.588011
## 281 282 283 284 285 286 287 288
## 5.809998 6.065841 5.670088 5.485997 1.720444 0.460659 0.004019 0.552030
## 289 290 291 292 293 294 295 296
## 1.153563 0.408696 1.278576 1.232461 1.232461 5.333110 5.333110 1.232461
## 297 298 299 300 301 302 303 304
## 0.657841 0.657841 1.223907 1.273789 0.903462 0.431246 0.388596 0.130004
## 305 306 307 308 309 310 311 312
## 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681
## 313 314 315 316 317 318 319 320
## 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681 0.730681
## 321 322 323 324 325 326 327 328
## 0.730681 0.730681 2.766052 4.215664 4.204283 4.204283 4.204283 5.022112
## 329 330 331 332 333 334 335 336
## 5.237510 5.237510 5.164118 2.788452 5.172184 4.838514 4.838514 4.838514
## 337 338 339 340 341 342 343 344
## 4.838514 2.039981 3.278505 5.533461 5.533461 5.533461 4.218526 4.218526
## 345 346 347 348 349 350 351 352
## 8.313752 8.313752 8.313752 8.642456 8.642456 8.642456 8.642456 8.868480
## 353 354 355 356 357 358 359 360
## 8.729123 8.015217 8.015217 8.015217 8.015217 8.015217 8.060607 8.015217
## 361 362 363 364 365 366 367 368
## 3.120847 3.014030 2.340587 2.257551 2.382146 0.661596 4.633715 2.110811
## 369 370 371 372 373 374 375 376
## 1.859324 1.877740 1.176540 10.823588 9.306410 7.854842 9.850179 10.769445
## 377 378 379 380 381 382 383 384
## 11.854208 11.578107 13.301104 13.850234 14.334685 10.815390 10.465704 10.017208
## 385 386 387 388 389 390 391 392
## 9.675993 9.353371 9.030021 13.460333 11.098649 10.827129 10.900202 10.571104
## 393 394 395 396 397 398 399 400
## 10.467844 10.384180 9.573346 8.310444 7.454699 9.414721 8.811303 9.272529
## 401 402 403 404 405 406 407 408
## 19.513793 20.272830 19.450144 20.560878 23.253747 28.368872 30.395795 25.079676
## 409 410 411 412 413 414 415
## 27.481017 41.577478 29.372089 30.523206 32.152446 31.107990 34.365793
Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.
We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).